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ABSTRACT 

We report a series of simulations of the formation of a star cluster similar to the Orion Nebula 
Cluster (ONC), including both radiative transfer and protostellar outflows, and starting from both 
smooth and self-consistent ly turbulent initial conditions. Each simulation forms > 150 stars and brown 
dwarfs, yielding a stellar mass distribution that ranges from < 0.1 Mq to > 10 Mq. We show that a 
simulation that begins with self-consistently turbulent density and velocity flelds embedded in a larger 
turbulent volume, and that includes protostellar outflows, produces an initial mass function (IMF) 
that is consistent both with that of the ONC and the Galactic fleld, at least within the statistical power 
provided by the number of stars formed in our simulations. This is the flrst simulation published to 
date that reproduces the observed IMF in a cluster large enough to contain massive stars, and where 
the peak of the mass function is determined by a fully self-consistent calculation of gas thermodynamics 
rather than a hand-imposed equation of state. This simulation also produces a star formation rate 
that, while still somewhat too high, is much closer to observed values than if we omit either the 
larger turbulent volume or the outflows. Moreover, we show that the combination of outflows, self- 
consistently turbulent initial conditions, and turbulence continually fed by motions on scales larger 
than that of the protocluster yields an IMF that is in agreement with observations and invariant 
with time, resolving the "overheating" problem in which simulations without these features have an 
IMF peak that shifts to progressively higher masses over time as more and more of the gas is heated, 
inconsistent with the observed invariance of the IMF. The simulation that matches the observed IMF 
also qualitatively reproduces the observed trend of stellar multiplicity strongly increasing with mass. 
We show that this simulation produces massive stars from distinct massive cores whose properties 
are consistent with those of observed massive cores. However, the stars formed in these cores also 
undergo dynamical interactions as they accrete that naturally produce Trapezium-like hierarchical 
multiple systems of massive stars. 

Subject headings: ISM: clouds — radiative transfer — stars: formation — stars: luminosity function, 
mass function — turbulence 



1. INTRODUCTION 

The origin of the the stellar initial mass function (IMF) 
is a classic problem in astrophysics. Since the IMF is 
most easily measured in young star clusters, and appears 
to be essentially the same in such clusters and in the fleld 
(e.g. Bastian et al. 2010), this problem is closely linked 
to the problem of how star clusters form. There have 
been numerous theoreti cal attacks on these twin prob- 
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lems (see the review by McKee fc Ostriker||2Q07 ), but a 
major breakthrough in the past lew years nas been the 
realization that the answer is tightly linked to the ques- 
tion of gas thermodynamics. An isothermal gas, even a 



or competitive accretion (e.g. Bonnell et al.|2UUla|bp , can 
predict a shape for the IMF, but in order to determine 
a characteristic scale must either appeal to additional 
physics or must deflne a flducial "cloud", whose mean 
density or other properties (e.g. the normalization of its 
linewidth-size relation) then determines the location of 
the IMF peak. In the latter approach, however, it is 
not clear on what scale one should measure cloud prop- 
erties: an entire CMC, with a mean density n ~ 10^ 



cm 



obeying the |Larson| ([1981| linewidth-size relation, 

10^ 



a massive clump with a mea n densit y n 



magnetized one, has no char acteristic mass scale (McKee 
et al.|2010[|Krumholz|2'QTT] ). This implies that the proE^ 
lem ol the origin ol the IMF, which is observed to be 
invariant in both its shape and its characteristic scale, is 
a separable one. 

Models that describe the behavior of a isothermal gas, 
such as those based on turbulent fragmentation (e.g. 
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et al.| |2003), or some other scaleY Ditterent choices yield 
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and a linew idth far above the Larson| value (e.g. [Shirley 



wildly varying characteristic masses. Moreover, cloud 
properties vary in sufficiently extreme galactic environ- 
ments, for example showing different linewidth-size rela- 
tions (e.g. Rosolowsky fc Blitz||2005 ). Despite this vari- 
ation, however, there is no evidence for a corresponding 
variation in the IMF. These problems strongly suggest 
that the origin of the IMF peak cannot be found in the 
physics of isothermal gas. Instead, models that seek to 
explain any characteristic mass scale in the IMF must 
appeal to departures from isothermality ("Rees'lQTGl 'Low! 
fc Lynden-Beh 1976, .Spaans fc Silk 2000; Larson„2005, 
|Krumholz||2011| ). 



Simulations of star formation mirror these trends de- 
pending on the physics they include. Isothermal simu- 
lations always produce a characteristic stellar mass that 
is determined b y the initial c onditions or the numerical 
resolution (e.g. Martel et al. 2006), and can always be 
rescaled to produce an arbitrary stellar mass scale. In 
contrast, those that include non-isothermality produce 
characteristic mass scales that are determined by the 
mechanism that causes them to depart from isothermal- 
ity, whether it b e an impos ed equation of state (e.g. Bate 

Bonnell 2005 Jappsen et al. 2005) or the inclusion of 
radiat ive transfer, either without (Bate 2009b, |2012[) or 



I 



with (Krumholz et al. 



2009 
stellar 



Urban et al. 
radiation. 



2007a 



|2010, ,2011, .Ottner et al.. 
20iop the further step of including 



Since comparisons between approxi- 
mate equations of state and radiative transfer calcula- 
tions show that the former offer only an extremely poor 
approximation, progress toward an understanding of the 
IMF's characteristic peak therefore requires radiation- 
hydrodynamic simulations. 

The radiation-hydrodynamic simulations that have 
been conducted thus far have demonstrated several 
promising features. First, radiation feedback suppresses 
the formation of brown dwarfs, reproducing the observed 
turn-down in the IMF at low masses (Bate 2009b Offner 



et al.||2009| ). Second, simulations includmg radiation 
feedback are able to suppress fragmentation in very dense 
regions, allowing the formation of massive stars when the 
conditions are approach those seen inrealregions of mas- 



sive star formation (Kr umholz et al.|2007a 2010). Third, 
radiative simulations produce an Iivll^' that does not vary 
with the properties of the star-forming cl oud in low mass, 
low-density environments (Bate 2009b| [2012 ), nor with 
the gas metallicity ( Myers et al.| |2011). 

While these results are encouraging, these simulation 
efforts have for the most part been limited either to sin- 
gle massive cores, or to clouds of l ow density and/or low 
mass. For example. Bate (2012) simulates a cloud of 
mean volume density 3 x 10"^ cm~^ and column density 
0.2 g cm~^ , forming ^ SO Mc?) o f stars, none larger than 
~ 3 Mq. Peters et al. (2010) do form massive stars, 
but from a cloud with a mean density of 10^ cm~^ and 
a column density of 0.026 g cm~^, far below the col- 
umn density at which radiative effects become impor- 



tant (Krumholz & McKee 2008 Krumholz et al. 2010) 



so low, in fact, as to be optically thin in the near- 
infrared. In contrast, the mea n mass and radius of the 
star- forming regions studied by Faiindez et al. (2004) im- 
plies a volume density > 10^ cm""^, a mass of ~ 5000 
M0, and a column density of 2 g cm~^, such that mul- 
tiple massive stars would be expected, and their radi- 
ation would be trapped effectively by the clou d's high 
optica l depth . Si milar Galaxy-w i de su rveys by |Shirley| 



et al] ([20031) and |Fontani et al.| (|2005|) that target re 
ns of 



gions 
erties. 
dN/dM (X M 



active star formation produce comparable prop- 
Indeed, the observed cluster mass funct ion is 
Lada||2003j [Fah et al. 



('Lada 



0. , ,, 

Chandar et al.pOfU), implying that a majority o: 



2009 



stars 
form in clusters larger than 1000 Mq in mass, large 
enough to possess O stars. The ONC, therefore, is a 
far more typical star-forming environment than most of 
the regions explored with radiation-hydrodynamic simu- 
lations thus far. 



In Krumholz et al. (2011 hereafter Paper I) we re- 
port eoTEeTTrsTTaHiaiioirTtydro dynamic simulations to 
probe this more typical regime of star formation; that 
calculation followed the collapse of a 1000 Mq cloud with 
a column density of 1 g cm~^, leading to the production 
of > 500 Mq worth of stars, with an IMF extending from 
~ 0.05 Mq brown dwarfs to ~ 30 Mq O stars. This 
calculation identified a problem. Radiative suppression 
of fragmentation, which seems necessary to explain the 
invariant peak in the IMF and avoid overproduction of 
brown dwarfs, became too efficient. As the calculation 
proceeded, the cloud underwent a global collapse, lead- 
ing to extremely high star formation rates and accretion 
luminosities. As a result, the gas heated up to the point 
where further star formation was suppressed. The net 
result was an IMF that was not invariant, but instead 
had a peak that moved to systematically higher masses 
as the calculation proceeded. At early times there were 
too few massive stars, and at late times too many. Since 
there is no plausible mechanism to guarantee that all 
star-forming clouds would stop producing stars at the 
same point in this evolution, this result was inconsistent 
with the observed universality of the IMF. 

In Paper I, we conjectured that the problem could be 
resolved by lowering the star formation rate per free-fall 
time, which would in turn lower the accretion luminos- 
ity. Such a change is required by observations even in 
the absence of the problems rapid star formation creates 
in the IMF, because observed star formation rates per 
free-fall time are always a few percent across a very wide 



range of star- forming environments (Krumholz fc Tan 
2007[ [Evans et al.||2009] [Krumholz et al.||2012[ ). In this 



paper we test that conjecture by performing additional 
simulations of the formation of ONC-like star clusters, 
with two extra pieces of physics that should lower the 
star formation rate per free-fall time. First, rather than 
simulating an isolated star-forming clump as in Paper I, 
we embed our initial clump in a larger volume of tur- 
bulent gas, and we initialize the simulations such that 
our clump has self-consistently generated turbulent den- 
sity and velocity structure. Second, we include proto- 
stellar outflows. A number of authors have shown that 
such outflows can inject significant energy into a star- 
forming cloud, driv ing its turbulence and lowering its 
star formation rat e (|Li fc N akamura 2006; Nakamur a fc 
Li 2007; Matzne?'2 007| | Wan g et al., |2010| [Cunningham 
et al.__2011J . Ideally, a third piece of physics should also 



be included: magnetic fields. Thes e both lower the sta r 



formation rate by themselves (e.g. [Price fc Bate 2009) 
and also enhance th e effectiveness of protostellar outflows 
(IWang et al.[ [2010 ). We plan to do so in future work. 

The remainder of this paper is organized as follows. In 
Section [2] we describe our numerical methods and simula- 
tion setup. In Section [3] we report the simulation results, 
and finally we discuss their implications and draw con- 
clusions in Section [U 



2. NUMERICAL METHODS 

We simulate star- forming clouds using the ORION code, 
which includes radiative transfer (jHowell fc Greenoughl 



2003 Krumholz et al. 



2QQ7bl [Shestakov fc^ ^ner 20 Q8J 



hydrodynamics ([Klein 1999p, self-gravity (['I'ruelove et al. 



1998[ ), accreting sink particles ( [Krumholz et al.[[2004[ ), 
and a model for protostellar evolution and feedback, in- 



eluding stellar radiation and outflows ( Offner et al.|2QQ9^ 
Cunningham et al.|2Qll ). Here we briefly summarize the 
equations we solve" the code itself, and the initial condi- 
tions for the simulations. For the first two of these topics, 
we refer the reader to Paper I for more details, since the 
physics included and the numerical methods are identical 
except where specified below. 

2.1. Equations and Algorithms 

ORION solves the equations of gravito-radiation- 
hydrodynamics in the two-temperature, mixed- frame 
flux-limit ed diffusion appro ximation. These equations 
are ( ,Krumholz et al.||2QQ7bl ) 
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In these equations, p, v, P, and e are the density, velocity, 
pressure, and total (thermal plus kinetic) energy density 
of the gas, E is the energy density of the radiation, (j) is 
the gravitational potential, kqy* and kqyi are the Planck 
and Rosseland mean opacities of the dust-plus-gas fluid, 
A is the flux-limiter, A is the rate of non-dust cooling 
(via line and continuum processes in gas at temperatures 

^ 10^ K where the dust sublimes), and m^ is the proton 
mass. For more information on the flux-limiter, hot gas 
cooling rate, and choice of dust opacities, we refer the 
reader to Paper I. 

Terms subscripted by i refer to stars; Xi^ Mi, and p^ 
are the position, mass, and momentum of the ith star. 



dt ' ' 






dt ' Mi 




-Pi = -MiV<P + Pi 




V2(/) = 47rG 


/9 + ^ Mi5ix - 


-X,) 




i 





Mi, Pi, and £i are the rate at which those stars add or 
remove mass, momentum, and energy from the gas, Li 
is the luminosity of star i, and Wi is the weighting ker- 
nel that spreads the stellar interaction over some number 
of computational cells. The equations we solve here dif- 
fer from those in Paper I in that, in addition to accretion 
(the terms subscripted with a), we also include protostel- 
lar winds (the terms subscripted with w). Stars accrete 
gas from the computati onal grid f ollowi ng the sink par- 
ticle method of Krumholz et aLj ( 2QQ4|) , and each sink 
particle is linked to a protostellar evolution code that 
computes the instantaneous stellar radius and luminos- 
ity based on the star's accretion hi story, follow i ng th e 
method described in the Appendix of Offner et al. ( 2009 ). 
In addition, during each time step, each star returns a 
portion of the mass it accretes to the grid in the form of 
a collimated protostellar wind. For details o f the numer- 
ical implementation, see Cunningham et al. (2011). Our 



wind parameters are the same as in that paper, i.e. each 
star ejects a fraction /^/(l + /w) = 0.21 of the gas it 
accretes (so /^ = 0.27), this material is launched with 
a velocity fy = 1/3 that of the Keplerian speed at the 
stellar surface, and the wind gas has a temperature 10^ 
K at launch. It is collimated along the axis defined by 
the stellar angular momentum vector. 

Orion solves Equations (fTl) - (|9]) within an overall 
adaptive mesh refinement (AMR) structure, in which the 
entire domain is discretized onto a coarse grid of size 
A^o cells on a side, denoted level 0. Sub-regions within 
the domain are then covered by progressively finer grids. 
The grid on level £ has a resolution a factor of 2^ better 
than that of the coarse grid, and evolves with a time 
step a factor of 2^ smaller. These grids are automatically 
added and removed on the fly as the calculation proceeds, 
based on user-specified criteria, up to some pre-specified 
maximum level L. 

2.2. Simulation Setup 

We compare three different simulations, which we re- 
fer to as smooth, no wind (SmNW), turbulent, no wind 
(TuNW), and turbulent, with winds (TuW); in terms of 
physics these differ in that the NW simulations have pro- 
tostellar outflows disabled. We summarize this and other 
properties of the simulations in Table [l] All simulations 
consist of a mass Mc = 1000 Mq of gas with a mean 
surface density He = 1 g cm~^ (arranged as described 
below). Throughout the cloud we set the gas tempera- 
ture and the radiation energy density to T^ = 10 K and 
E = aTg = 7.56 x 10~^^ erg cm~^, respectively. 

For run SmNW, we use a setup identical to run HR 
from Paper I (though here we have continued the simu- 
lation further in time than we described in that paper), 
so we only briefly discuss its properties here, and re- 
fer readers to Paper I for a fuller description. In run 
SmNW, the initial gas distribution is a sphere with a ra- 
dius Re = 0.26 pc. The density distribution is smooth, 
and consists of a central core of uniform density that ex- 
tends to half the cloud's radius, surrounded by an outer 
region within which the density falls off with radius as 
as suggested by observations of massiv e clumps 



-1.5 



(e.g. |Sridharan et al.|200"5l|Beuther et al.l2006D . The gas 
is given an initial turbulent velocity tield with a disper- 
sion of (Jc = 2.9 km s~^ (one-dimensional), correspond- 
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No 


1000 
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Fig. 1. — Column density distribution in the turbulent initial 
conditions used for runs TuNW and TuW. 

ing to an initial virial ratio a = bcr^Rc/GMc = 2.5. The 
velocity power spectrum is P{k) oc /c~^, drawn without 
imposing any bias in favor of sol enoidal or compressiv e 



modes following the procedure of Dubinski et al 



( 2Q12| )p] 



Note. — Col. 3: cloud mass. Col. 4: cloud radius (for run SmNW) or box size (for runs TuNW 
and TuW). Col. 6: mass weighted-mean density at time t = 0. Col. 7: free-fall time computed 
using {p)m- Col. 8: size of computational box. Col. 9: number of cells per linear dimension on 
the coarsest AMR level. Col. 10: finest AMR level. Col. 11: grid resolution on the finest AMR 
level. 

In runs TuNW and TuW we initialize so that, unlike 
in run SmNW, both the initial density and velocity fields 
are self-consistently turbulent. We set up a periodic do- 
main of length ic = 0.46 pc on a side, so that Sc = 1 
g cm~^ averaged over the box. To initialize the simula- 
tion, we impose the same turbulent velocity field as in 
run SmNW, scaled to a velocity dispersion ac = 1.4 km 
s~^, corresponding to a = 1/2 if we use -^c/2 in place of 
Re. Although this means the gas is less turbulent ini- 
tially than in run SmNW, as we see below, damping of 
the turbulence in run SmNW brings the a values closer 
together as the runs progress. To produce a density field 
consist with this velocity field, we drive the turbulence 
and allow the simulation to evolve for two crossing times. 
During this period we turn off both gravity and radiation, 
and we hold the gas isothermal at a temperature T^ = 10 
K by setting the gas ratio of specific heats to 7 = 1.0001; 
since, in the absence of stellar sources, molecular cloud 
gas is close to isothermal, this should be a very good 
approximation, and ignoring radiation during this setup 
phase significantly reduces the computational cost. Dur- 
ing this setup phase we also fix the computational reso- 
lution at 512^ cells, with no further refinement. At the 
end of two crossing times we turn off driving, change the 
gas ratio of specific heats to 7 = 5/3, turn on gravity and 
radiation, and return to our normal refinement criteria 
(see below). This state represents the initial condition 
for runs TuNW and TuW. Note that, since the turbu- 
lence is driven mostly on large scales, the result of this 
procedure is essentially a single, dense, turbulent cloud, 
surrounded by lower density turbulent material; we show 
this state in Figure [T] This clump is therefore analogous 
to the isolated one in run SmNW, but is surrounded by a 
realistic turbulent environment rather than an artificial 
hot ambient medium. 

In all simulations the refinement criteria used to add 
higher resolution grids are the same. Specifically, we 
add resolution in any cell that satisfies one of the fol- 
lowing three conditions: (1) the density in the cell ex - 
ceeds the local Jeans density (Truelove et al. 1997), 
Pj = J'^TTc'l/GAxf, where J = 1/4 is the Jeans number, 

Cs = ^s/kBT/jii is the isothermal sound speed, and Axi is 
the grid spacing on AMR level /; (2) the radiation energy 
gradient is sharp enough so that \\/E\/E > 0.15/ Axi (al- 
though we sometimes temporarily reduce the coefficient 
below 0.15 for stability reasons in the TuNW and TuW 
runs); (3) the cell is within a distance of WAxi of any 
star particle. We refine to a maximum resolution of 49 
AU (L = 5) in run SmNW, and 23 AU (L = 4) in runs 



(1995) 



Outside the sphere of gas we place a zero-opacity am- 
bient medium with a temperature 100 times larger and 
a density 100 times smaller than that of the gas at the 
sphere's edge. We emphasize that, because the density 
gradient in the gas only extends to half the initial ra- 
dius, the overall center to edge density contrast is only 
a factor of 2.8, substantially less than that induced by 
the turbulent shocks. Thus this initial condition is quite 
similar to that adopted b y other authors who h ave simu- 
lated isolated clouds 



e.g. 



Bonneh et al. (2003) and Bate 



^ An additiona l diff erence betwe en our setup and that of |Bon-| 
|nell et aT] (2003) and Bate| ( |2012| ) is that we place an ambient 
medium outside our cloudmat is m thermal pressure balance with 
the material at the cloud edge, while the smoothed particle hy- 
drodynamics (SPH) simulations of Bonnell et al. and Bate have 
a vacuum outside their clouds. However, this ditterence is almost 
certainly negligible. The thermal pressure of our ambient medium 
is set equal to the thermal pressure of the cloud, which is smaller 
than either the ram pressure or the self-gravitational weight of the 
cloud by a factor of ~ 100. Thus the extra pressure provided by the 
external medium will enhance the collapse that would occur due to 
gravity alone by only ~ 1%. Even this is likely an overestimate of 
the difference between the two simulation methods, because, while 
formally the SPH simulations have vacuum outside their clouds, 
SPH creates an artificial surface tension at density discontinuities 
(fPrice 2008), and this will act very much like a confining external 
pressure. Our Eulerian simulation method does not suffer from 
this problem. 



TuNW and TuW. Finally, we note that, while the hy- 
drodynamic and gravitational boundary conditions are 
necessarily different in the smooth and turbulent runs, 
the great majority of the star formation in the turbulent 
runs occurs in subregions much smaller than the entire 
computational volume, and thus the periodic boundary 
conditions have minimal impact. We also impose Mar- 
shak boundary conditions on the radiation in runs TuNW 
and TuW in order to let radiation es cape the computa- 
tional volume (ci. [Offner et al.||2QQ9D . 



3. RESULTS 

Before examining the results of our simulations, we 
first mention two subtleties in the analysis that apply 
to the remainder of this discussion. First, since we are 
comparing runs with different initial conditions, it is im- 
portant to normalize the times so that differences be- 
tween the runs reflect the underlying physical behavior, 
and not simply that the dynamical time is different in 
different cases. Moreover, in the runs with turbulent 
initial conditions, the strong initial turbulence guaran- 
tees that the majority of the mass is compressed into 
structures that are significantly denser than the volume- 
averaged density. Given these considerations, the most 
natural approach, which we adopt, is to measure times 
in units of the free-fall time tff = Y^37r/32Gp evaluated 
at a density equal to the initial mass-weighted density 
{p)mi since this is the dynamical time appropriate to the 
bulk of the matter. This approach also has the advan- 
tage that it is the most natural basis for observational 
comparison, since an observation would detect the bulk 
of the mass, and would be sensitive to the typical den- 
sity at which this mass resides. For this reason, in what 
follows whenever we refer to times, we normalize to tff 
defi ned in this manner. We report this quantity in Table 



3 



The second subtlety is that, as in Paper I, we only 
regard stars as collapsed objects once their mass exceeds 
0.05 Mq^ based on one-dimensional calculations of the 
mass at which second collapse to stellar densities occurs 



(Masunaga & Inutsuka 



2000). We allow smaller objects 



to merge with one another and with more massive stars. 
We therefore restrict our analysis to objects larger than 
this mass. 



3.1. Overall Evolution and Morphology 

In Figures [2J [3J and [4] we show the large-scale column 
density and density-weighted temperature distributions 
in for runs SmNW, TuNW, and TuW, respectively. In 
all three we observe the same general trend: the turbu- 
lence creates an overdense region, which then begins to 
collapse and form stars. The collapsing structures are fil- 
amentary, and the stars are born along the filaments, and 
particularly at the nodes where the filaments intersect. 
The temperature is initially small, but as stars form hot 
spots around individual stars appear, and these gradu- 
ally spread over timejj 

^ Note that in Paper I we instead used the volume weighted 
mean density to compute the free-fall time for run SmNW; how- 
ever, because the initial density field is very smooth, the difference 
between volume- and mass-weighted mean density free-fall times 
for this run is only ~ 20%. 

^ In runs TuNW and TuW there are sometimes brief increases in 
the overall background temperature level visible in some snapshots. 






T 



Fig. 2. — Column density (left) and density- weighted mean tem- 
perature (right) in run SmNW. The times of each pair of images 
are indicated in the right column, running from t/% = to 1.25 
in steps of 0.25. In the column density plot, white circles indicate 
the positions of star particles, with the size of the circle indicating 
the mass of the star. In the right column, the temperature shown 
is the radiation temperature Tr^ defined implicitly hj E = aT^ . 
We show this rather than the gas temperature because the gas and 
radiation temperatures are nearly equal everywhere in the cloud, 
except in the hot ambient medium outside the cloud in run SmNW, 
and in material ejected by protostellar outflows in run TuW. Us- 
ing the radiation temperature provides a convenient means to fllter 
this contribution. 

There are a few interesting points to take from these 
plots. One involves the morphology of the heated regions. 
In the turbulent runs, even at late times the temperature 
distribution looks more like a series of islands of heated 
gas surrounded by a large medium that is either at or 
quite near to the background temperature. In effect, one 
can discern something like individual protostellar cores 
that are heated by the star or star system embedded 

but these are short-lived and small, generally keeping the temper- 
ature < 20 K. These flashes are associated with brief increases in 
the accretion luminosity that are large enough to heat the entire 
simulation volume above 10 K for short periods. 
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Fig. 3. Same as Figure [2] but for run TuNW. Note that 

the color scales are the same, out the size of the region shown is 
slightly different. 

within them. In contrast, by the end of run SmNW there 
is simply a single, concentrated region of heating, and 
one cannot discern individual cores any more. As we 
show below, this difference proves to be important in 
determining the evolution of the IMF. 

Another interesting point is that the overall morphol- 
ogy is surprisingly similar in runs TuW and TuNW, de- 
spite the change in whether we include protostellar winds 
or not. Partly this is a function of the fact that wind- 
blown bubbles are fairly low column-density structures, 
and that we are looking at static slices. In an anima- 
tion of the column density field, one readily discern out- 
flows driving shells of gas orthogonal to the filaments. 
However, this clearly has a relatively small effect on the 
large-scale morphology. 

3.2. Star Formation Rate and History 

In Figure [5] we show the star formation history of each 
of our simulations. The most immediate and striking 
thing about the Figure is the difference in star forma- 




FlG. 4. — Same as Figure |3] but for run TuW. 




Fig. 5. — Total mass in stars (top) and total number of stars 
(bottom) as a function of time in runs SmNW, TuNW, and TuW. 





TABLE 2 

Simulation Outcomes 




Name 


tfin/tff M,^^^/Mc 


A^* 


lO^M* 

(Mo yr-i) 


eff 


SmNW 

TuNW 

TuW 


1.25 0.70 
1.39 0.20 
1.32 0.15 


540 
127 

158 


16 

7.4 
6.2 


1.78 
0.33 
0.28 



Note. — Col. 2: time at which run was stopped. Col. 3: 
total stellar mass at the end of the run. Col. 4: number of 
stars present at the end of the run. Col. 5: time-averaged 
star formation rate in the run, measuring from formation 
of the first star to the end of the run. Col. 6: dimensionless 
star formation rate e^ = M^ /[(I / 2)Mc / 1^]. 

tion histories between the smooth and turbulent runs. 
Run SmNW starts off its star formation more slowly 
than TuNW or TuW, which is not surprising since its 
initial density structure is smooth, and possesses no high 
density peaks that collapse quickly. However, star for- 
mation in that run accelerates dramatically as the time 
approaches tff. Late in the simulation, the star forma- 
tion rate approaches ^ 1 — 2Mc/tf[. In contrast, in runs 
TuNW and TuW the star formation rate is roughly con- 
stant and fairly low. After a time tff, only about 10% of 
the mass has been turned into stars. There is no obvious 
acceleration with time. Note that, although part of the 
difference in star formation rates comes from the differ- 
ence in free-fall times between the turbulent and smooth 
runs, even if we were to measure in seconds rather than 
free-fall times run SmNW would have a much larger star 
formation rate than TuNW or TuW. We summarize the 
dimensional and dimensionless star formation rates in 
the simulations in Table [2] Observationally, the dimen- 
sionless star formation rate (Krumholz fc McKee,2QQ5J 



eff 



M* 



(l/2)M,/tff' 



(10) 



where the factor of 1/2 arises because half the cloud mass 
is above the density {p)m used to define tff , is ~ 1%, with 
roughly half a dex scatter, over a very wide ran^e of den- 
sities and galactic environments ferumholz 



Tan2007 



Evans et al. 



2009 



Krumholz et al. 



2012 



I" J Comparing 
W and TuW is 



to the table, we see that Cff in runs I'uN 
still roughly an order of magnitude too high compared 
to observations, but it is roughly an order of magnitude 
lower than in run SmNW. We discuss the origin of the re- 
maining discrepa ncy b etween TuW and the observations 
further in Section lOl 

The difference in star formation rate (SFR) between 
the runs may be understood readily if we consider what 
happens to the turbulence, which, in these runs with no 
magnetic fields, is the main mechanism for regulating 

^ There is some subtlety in the observational comparison here, 
because real observations usually have an upper limit on the den- 
sity to which they are sensitive, for example because the tracer 
being used depletes or becomes very optically thick at high den- 
sity. Since we include all the mass above {p)m in our computation 
of eff , we are not capturing this effect. However, the change in mass 
it would induce is small, because both real star-forming clouds and 
our simulated turbulent clouds have density probability distribu- 
tion functions that are sharply declining at densities above the 
peak. Thus the amount of mass missed due to the density upper 
limit in the observations is likely to be very small. 



the SFR. In run SmNW, the initial turbulence present 
in the gas decays, and after one crossing time, which is 
~ tff, this decay gives rise to a global collapse and an 
accelerating star formation rate. In contrast, for runs 
TuNW and TuW, the box crossing time, and thus the 
turbulent decay time, is significantly longer than the 
free-fall time at the mass-weighted mean density. It is 
comparable to the free-fall time at the volume-weighted 
mean density, which is much longer. The difference in 
star formation history between the runs makes a criti- 
cal point: it matters for their star formation rates that 
star-forming dense clumps like the one out of which the 
ONC formed are not isolated objects. They are instead 
the inner parts of larger turbulent structures, and the 
energy from those larger scales is able to cascade down 
to smaller scales and maintain the turbulence for longer 
than the dynamical times of the small clumps. The tur- 
bulent decay timescale in a proto-ONC gas clump is the 
crossing time of its parent molecular cloud, not the cross- 
ing time of the small clump. This point has previou sly 
been made by iFalceta-Goncalves & Lazarian (2011) in 
the context of non-sell-gravitating turbulence, and our 
work strongly confirms their conclusion and extends it 
to the self-gravitating case. 

In contrast, the differences between the two turbulent 
runs are relatively small. The star formation rate mea- 
sured by mass (as opposed to number of stars) is ~ 20% 
lower in TuW than in TuNW. Since our model for pro- 
tostellar outflows prescribes that 27% of the mass that 
reaches a star particle (and thus the inner wind launching 
region) be ejected in an outflow, this reduction in the star 
formation rate is surprisingly small. This implies that 
there must be very little entrainment of additional mate- 
rial by the outflows, and even that some of the material 
that is entrained by outflows must be rapidly stopped 
and recycled back into the star-forming region. Visual 
inspection of the morphology of the outflows and accre- 
tion flows confirms that this is in fact the case: outflow 
shocks visible in the animations are generally traveling 
at right angles to the filaments feeding the stars. This is 
not an accident. Each star launches its bipolar outflow 
along the axis specified by its angular momentum vec- 
tor. If stars are being fed primarily by filaments lying 
in a plane, as is the case in all our runs, then most of 
their angular momentum vectors tend to be perpendicu- 
lar to that plane, producing relatively little entrainment. 
The minority of outflows that do end up aligning with 
the filaments possess too little momentum to significantly 
hinder the accretion flow, and the matter they do eject is 
stopped by the greater ram pressure of the infalling gas. 
As a result, it is re- accreted fairly rapidly. Whether this 
behavior is actually realistic is a separate question, one 
to which we return in Section [4j 

3.3. The IMF 

Another interesting feature in Figure [5] is that the mass 
in stars in run TuW is ^ 20% smaller than in TuNW 
at equal times, but that the number of stars is ~ 20% 
larger in TuW. This indicates an important shift in the 
stellar IMF between the runs. Although it is less obvious 
visually from Figure [5| there are also very important dif- 
ferences in the IMF between run SmNW and the other 
two runs. We now examine these. 

For the purposes of quantitative comparison between 
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Fig. 6. — Evolution of the IMF over time in the three simula- 
tions. Thick lines indicate the 50% percentile mass M50 (see main 
text for formal definition), while the shaded regions indicate the 
range between the 25th and 75th percentile masses M25 and M75. 
Thin lines indicate the mean mass M. Colors indicate the run, as 
described in the legend. Circles long the thick lines indicate the 
points at which the stellar mass reaches 50 Mq, 100 Mq, 150 Mq, 
etc. For comparison, thick gray unbroken horizontal lines show 
^25, M50, and M75 for a fully sampled Da Rio et al. (2012) IMF 
(Equation |13| ), and the thin gray unbroken horizontal line shows 
M for this IMF. Dashed lines show the equivalent quantities for a 
|Chabrier| ( [2005| ) IMF. 

different simulations, and between simulations and ob- 
servations, it is helpful to examine percentiles in the cu- 
mulative mass distribution function for stars produced in 
the runs. We define the nth percentile mass M^ implic- 
itly via the equation 



E 



100 



E^ 



(11) 



where m^ is the mass of each individual star, the first 
sum runs over stars with masses m^ < Mn^ and the 
second sum rus over all stars. Thus, for example, M50 
is defined by the condition that sum of the masses of 
all stars smaller than M50 constitutes 50% of the total 
stellar mass. We also examine the mean stellar mass, 
defined by 

- £^ (1.) 



M 



where N^ is the total number of stars. We can measure 
each of these quantities directly from our simulations at 
every time. We can also compare the simulation IMFs 
to observed ones. We select two observational IMFs for 



comparison. In ONC, Da Rio et al. (2012) find for low 
mass stars an IMF well-tit by a lognormal function with a 
width of <j = 0.44 in log m^ , centered on log m^^c = —0.45 
( measured in Mq ; their Table 3). The highest mass bin 
in |Da Rio et al. s sample is ^ 2 M(:) , so to extend this 
to higher masses we adopt a [Chabrier (2003) functional 
form in which the lognormal at low mass has a powerlaw 
tail of slope —1.35 at high mass. Thus the observed ONC 
IMF to which we compare is 



dN 



d\ogm^ 



oc 



-(logm,-logm*,e)V2a'^ m^ < M(7 

m* > M0, 



-logmi^^/2a^ 



-1.35 



(13) 

where all masses are in Solar units, over a range from 
0.05-150 Mq. For this IMF, M25 = 0.69 M©, M50 = 1.8 



comparison IMF is the system IMF of Chabrier (2003 



ai 



2005[) for the galactic field, which also seems to tit other 



star clusters reasonably well ([Parravano et al.|2Qll[ ) . We 
use the system rather than the single star liVll^' because 
we do not resolve tight binaries. This IMF has the same 
functional form as Equation (13), but with logm^^c = 
—0.60 and a = 0.55. The corresponding percentile and 
mean values are M25 = 0.63 M©, M50 = 1.7 M^t) , M75 = 
7.2 M0, and M = 0.73 Mq. The 'Chabrier' and JDa Rio 
|et al., IMFs differ significantly in the number of brown 
dwarfs and very low mass stars they predict, but converge 
at masses above a few tenths of Mq . For a discussion of 
possible origins of t he discrepancy betwe en the two IMFs, 
we refer readers to Da Rio et al. (20121). 

In Figu re [6| we plot the time evolution of M25, M50, 
M75, and M in each of our simulations, and for the ob- 



served |Da Rio et al] |2012D and |Chabrier| ( |2005p IMFs. 
The figure immediately reveals some mterestmg results. 
First, we see that in run TuW the IMF is in remarkably 
good agreement with the observed IMFs. At the end of 
the si mulations, the mean mass M agrees with the ob- 



served iDaRioetaLj value to better than 20%, and the 
50th percentile mass M50 to less than a factor of 2; we 
show below that this level of disagreement is consistent 
with coming simply from statistical sampling variance. 
Moreover, and perhaps more importantly, the agreement 
is good at almost all times when there is a significant 
mass of stars present, because the IMF in run TuW is 
very stable over time. From ~ 0.7tff, when the total 
stellar mass reaches 50 M©, to ~ 1.3tff, when it reaches 
150 M0, we find that M50, M25, and M stay constant 
to within ~ 50%; M75 changes slightly more, almost cer- 
tainly as a result of under-sampling the high end of the 
IMF when there are relatively few stars. The change be- 
comes even smaller at later times. From the time when 
10% of the mass is in stars {t ~ l.Otff ) to when 15% is in 
stars, M50 changes by less than 5% and M by less than 
10%. 

In contrast, for run TuNW the mean mass is relatively 
stable, but M50 rises systematically with time, increasing 
by a factor of 2.2 as the stellar mass grows from 50 to 200 
M0, corresponding to times t/tf^ ~ 0.7 to 1.4. This re- 
flects more rapid growth of the more massive stars in the 
run where winds do not suppress accretion. Moreover, in 
this run the rate at which new stars form is lower than in 
run TuW. The agreement with observations in this case 
is clearly weaker; the run produces an IMF that is too 
top-heavy. 

The changes with time in run TuNW, however, are 
small compared to those that occur in run SmNW. There 
M50 and M increase by nearly an order of magnitude in 
a time less than 0.5tff. Each increase in stellar mass of 
50 Mq is accompanied by a factor of ~ 2 gain in M50. 
This pattern of growth occurs due to the "overheating" 
problem discussed in Paper I: in run SmNW, star for- 
mation is much too rapid and too concentrated, and this 
produces a rapidly rising accretion luminosity that heats 
the gas mass to the point where the Bonnor-Ebert mass 
is too large for stars for new small stars to form. Accre- 
tion continues, but it is entirely captured by the existing 
stellar population, leading to an IMF whose mean and 
median mass rise with time. Moreover, since all the stars 
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Fig. 7. — Cumulative mass functions in the simulations, com- 
pared to observations. The top set of panels shows the cumulative 
distribution by mass, and the bottom shows the distribution by 
number. Within each set of panels, columns show the results from 
runs SmNW, TuNW, and TuW, as indicated. Rows correspond to 
times t/tff = 0.75, 1.0, and 1.25, as indicated. In each panel, the 
colored line indicates the fraction of stellar mass /m(< ^^*) or 
the fraction of the number of stars /Ar(< fn^) in stars with mass 
less than vn^ in the simulation, and the gray band indicates the 
range from the 10th to 90th percentile resulting from drawing a 
large number of clusters from the Da Rio et al. ( 2012j) IMF. The 
hatched band is the 10th to 90th percentile range tor a |Chabrier| 
f2005) IMF. For details on how this drawing is done, see the 
Appendix to Paper I. The label in each panel indicates the total 
mass or number of stars at that time in that simulation. 



are growing in lockstep the mass distribution in this run 
is too narrow as weh. 

We can also make this comparison more quantitatively. 
In Figures [7| and [8] we show the cumulative and differ- 
ential mass distributions produced in our simulations at 
various times, and compare to observed IMFs. At each 
time in the simulations, we can quantitatively described 
the level of consistency or inconsistency between the sim- 
ulated and observed IMFs using a Kolmogorov-Smirnov 
(KS) test. We plot the result in Figure [9J Examining the 
figures, we see that run SmNW is strongly inconsistent 
with both observed IMFs at most times. At early times 
the IMF is too bottom-heavy, but as time increases the 
IMF peak shifts to higher masses. Around tji^ = 1 run 
SmNW is fully consistent with the Chabrier IMF, and 
marginally consistent with the da Rio one, but at later 
times the IMF peak continues to shift to higher values 



Fig. 8. — Same as Figure |7| but showing differential rather than 
cumulative mass distributions. The histogram value in each bin 
shows the total fraction of all stellar mass (for the top panels) or 
the total fraction of the number of stars (for the bottom panels) 
falling within that bin. Thick colored lines indicate the simulation 
result, and gray lines indicate the results of drawing an equal 
mass stellar population from the Da Rio et al.| (I2012| IMF. For 
the gray histogram, the histogram values give the median result, 
and the vertical lines indic ate the range fr om the 10th to the 
90th percentile. We omit the |Chabrier| ( pOOSt IMF here to reduce 
clutter. 



and becomes inconsistent with both once more. This is 
the overheating problem described in Paper I. For run 
TuNW the IMF peak does not shift systematically with 
time, and so there is no overheating problem. However, 
the absolute value of the mean mass is systematically too 
high, as shown in Figure [6J As a result, the overall level 
of agreement between the simulation and the observed 
IMFs is poor. On the other hand, run TuW is generally 
statistically consistent with both the da Rio and Chabrier 
IMFs at most times. 

It is important to add some caveats to this result. 
First, the KS statistic does not account for the observa- 
tional and systematic uncertainties in the observed IMFs. 
Were these uncertainties to be included, it is entirely pos- 
sible that run TuNW would be consistent within them, 
and perhaps even that run SmNW would be, at least 
for a longer period of time. Second, the KS test itself is 
an imperfect tool. It is most sensitive to differences in 
distributions near the 50th percentile, and less sensitive 
to differences on the tail of the distribution. Thus, for 
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Fig. 9. — Level of statistical agreement between the simulation 
and observed IMFs as a function of time. At each time, the 
quantity plotted is the result of a KS test comparing the three 
simulations (green for Sm NW, b lue of for TuNW, and red for 
TuW) to the D a Rio et aI1 ( |2012| thick lines) and Chabrier (2005 
thin lines) IMi^'s^ i'Jie rignt axis shows the P-value returned by 
the KS test, where 1 — P is the confidence level at which we can 
rule out the null hypothesis that the simulation and observed 
IMFs are drawn from the same underlying distribution. The left 
axis shows the equivalent confidence level measured in number of 
standard deviations, which is related by N^ = \/2 eric~^ {P) ^ with 
erfc the complementary error function. The dashed horizontal 
black line indicates a confidence level of 2.5cr. 

example, Figure |8] shows that run SmNW has a shght 
excess of stars in the ~ 3 — 10 M© range that is clearly 
visible in a differential mass function on a logarithmic 
axis. The KS test does not regard this excess as sta- 
tistically significant, but it is conceivable that a more 
sensitive statistical test might. Indeed, one can get a 
sense of the level of statistical power that the KS test 
provides when applied to our simulations from the fact 
that, formally, our simulations are consistent with both 
the da Rio and Chabrier IMFs. This is partly because 
we are not performing a comparison in the mass range 
0.01 — 0.05 Mq where the two distributions are most dif- 
ferent, but it is also partly because, with only 158 stars 
in run TuW, there is significant sampling noise. 

3.4. Gas Thermodynamics 

The fragmentation of the gas is driven by its thermo- 
dynamics, and we can gain insight into the differences 
in outcome between the runs by examining the temper- 
ature structure of the gas. In Figure [lO] we show phase 
diagrams of the three runs at three diirerent times. Not 
surprisingly, each of the runs is quite different. First ex- 
amining run SmNW, we note that, at time t/tf^ = 0.75, 
the bulk of the gas in run SmNW is cooler than in the 
other two runs. This reflects the fact that the total stellar 
mass in run SmNW is comparable to that in runs TuNW 
and TuW at this point. Since the free-fall time is longer 
in run SmNW, this corresponds to a lower total accretion 
rate and thus a lower accretion luminosity. However, as 
star formation in run SmNW accelerates, the accretion 
luminosity rises and the gas heats, while the gas in the 
other two runs stays relatively cool. Quantitatively, at 
the final times shown in the bottom row of Figure 10 , 42% 
of all the gas is at temperatures above 50 K in run SmNW 
(excluding the ambient medium); the equivalent Figures 
in both runs TuNW and TuW are 7%. It is important to 
note that this difference is driven by accretion luminosity 



and not by the intrinsic luminosity of massive stars. If 
we instead examine run SmNW at time t/tf[ = 1.0, the 
most massive star present is 8.8 Mq, smaller than the 
most massive stars present at time t/tf^ = 1.25 in runs 
TuNW (13.3 Mq) and TuW (9.9 M©). Nonetheless, we 
still find that 23% of the mass is at temperatures above 
50 K, and 64% is above 30 K, i.e. there is more hot gas 
in run SmNW even when the individual stars are less 
massive. 

The rapid heating in run SmNW gives rise to the over- 
heating problem identified in Paper I - bulk heating of 
all the gas makes it impossible for small stars to form, 
thus shifting the IMF systematically to higher mass as 
time goes on. Runs TuNW and TuW clearly do not suffer 
from this problem. Even at late times, the great majority 
of their gas is at temperatures of no more than 10 — 15 K, 
and there is very little material at temperatures of more 
than 50 K. Although there clearly is gas being warmed by 
stars in these runs, there remain pockets of cold, gas at 
densities > 10~^^ g cm~^ and temperatures < 15 K that 



is capable of producing new stars with masses 
Mq. These are visible in Figure 



10 



0.01 
where the phase 
diagram reveals the presence of material for which the 
Bonnor-Ebert mass, 



Mbe = 1.18- 



(14) 



is below 0.01 Mq. In contrast, at late times in run 
SmNW there no material for which Mbe is this small. 
To be quantitative, at the times shown in the final panel 
of Figure 10, run SmNW contains only 1.8 x 10~^ Mq 
of materiaTm the density and temperature region where 
Mbe < 0.01 M©, i.e. too little mass to actually create 
a star. The corresponding figures for runs TuNW and 
TuW are 0.49 and 1.0 M©, respectively, making it pos- 
sible for new brown dwarfs to form. It is interesting to 
note that the amount of cold, high-density gas is gener- 
ally greater in run TuW than in run TuNW. This is likely 
an effect of the reduced accretion rate and changed IMF 
in run TuW compared to TuNW, both of which serve 
to generally lower the accretion luminosity and thus the 
heating rate. The spatial distribution of the star forma- 
tion may also play a role: in run SmNW, because there is 
no pre-existing density structure at the start of the sim- 
ulation and because the turbulence decays rapidly, all 
the stars and gas become concentrated in a single dom- 
inant cluster, where stellar heating is very intense. In 
runs TuNW and TuW, the combination of a pre-existing 
density structure present in the initial conditions and 
the non-decay of turbulence throughout the simulation 
serves to break star formation up into several subclus- 
ters, within each of which stellar heating is less intense. 

3.5. Massive Cores and Massive Stars 

Run TuW is the first published simulation that in- 
cludes radiative and protostellar outflow feedback, pro- 
duces an IMF that is in good agreement with the ob- 
served IMF over a broad mass range, and forms a large 
enough cluster for there to be massive stars present. It 
is therefore important to pay particular attention to the 
processes by which those massive stars form. We turn 
now to the properties of the massive stars in run TuW. 

There has been considerable discussion in the litera- 
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Fig. 10. — Phase diagrams of the three runs at different times. The three columns correspond to runs SmNW, TuNW, and TuW, as 
indicated. The three rows correspond to times t/% = 0.75, 1.0, and 1.25. In each panel, the color indicates the gas mass in a given bin of 
density and temperature; bins are 0.025 dex wide in both p and T. The color scale is normalized so that the bin containing the largest 
amount of mass is 1.0. The long-dashed line indicates the locus in density and temperature at which the code inserts sink particles. 
The short-dashed lines indicate the locus in density and temperature where the Bonnor-Ebert mass is 0.01 M©, 0.1 Mq, and 1 Mq as 
indicated. Note that gas in the winds is run TuW is heated to ~ lO'^ K, well above the temperature range shown here, but there is 
relatively little mass at these temperatures. 



ture about whether massive st ars form from distinct mas- 
sive protostehar cores ( Padoan 1995| iPado an fc Nordlund 
2002] IMcKee k Tan|r2 UU2 , 2003; Krumhoi ret al l|i^007al 
Hennebelle fc Teyssier; |2QQ8t [Hennebelle fc Ch abrier 



2009), or whether all stars are born from cores with 
masses < 1 M©, and massive stars subsequently grow 
from these smal l seeds by Bond i-Hoyle accretion (Bon^ 
neh et all [19971 |2001a|b[ [20041 ["2 006; Bonneh & Bate 
2002^ ,2006i [Bate fc Bonnell||2005[ [Smith et aL„2009a,b^ . 



A number ot authors have also proposed hybrid models , 
in which massive stars form from gravitationally bound 
gas structures, but these structures are assembled and 
fed from larger scales at the same time as they form 
massive stars ( .Peretto et al.|2006[[Wang etal.| 20101 ). To 
address this question, we examme the four most massive 
stars present at the end of run TuW; these have masses 
of 10.8, 9.8, 8.8, and 8.3 M©, respectively, and thus each 
is large enough that, even if it were to accrete no fur- 
ther, it would be expected to end its life as a supernova. 
For comparison, we also examine the four stars whose 
masses are closest to the median mass at the end of the 
simulation, 0.34 Mq. For each of these stars, we identify 



the point in space and time at which that star first ap- 
peared in our simulations, and examine the gas density 
distribution in its vicinity. 

We show the results in Figure [TT] for the high mass 
cores and Figure [T2| for the low mass cores. To facili- 
tate comparison with observations, in addition to show- 
ing the true gas density distribution, we show the dis- 
tribution smeared with a 1700 AU Gaussian beam; we 
choose this size scale because it is approximately the spa- 
tial resolution of the highest published resolution maps 



Bontemps 



of massive cores (e.g. Beuther & Schilke 2004 
et al. 2010), though the Atacama Large Millimeter Ar- 
ray (ALMA) will soon produce images at significantly 
higher resolution. Figure 11 demonstrates that the mas- 
sive stars in our simulation form in distinct, massive 
overdensities that can be identified as cores. Their char- 
acteristic sizes, determined from visual inspection, are 
roughly 0.01 pc. Comparing the gravitational and kinetic 
energies in this structures shows that they are roughly 
gravitationally bound and virialized. The fiows within 
them are highly supersonic, producing a filamentary mor- 
phology. Nonetheless, these objects are not highly sub- 
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Fig. 11. — Images of the initial cores that produced the four most massive stars in simulation TuW. In each column, the upper image 
shows the column density distribution, centered on the ~ 0.05 Mq protostars that will grow to be massive stars. The lower image shows 
the same column density distribution, smeared with a 1700 AU Gaussian beam. In the upper panels we indicate the mass of the core 
(defined as the projected mass within a radius of 0.01 pc, as indicated by the dashed circles) and the time at which the snapshot is taken. 
In the lower panels we indicate the core mass that would be inferred from the beam-smeared image and the final mass of the resulting star. 
Note that the second and fourth columns are nearly identical because two of the final massive stars both form in the same core. 
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Fig. 12. — Same as Figure [TT] but for the four stars closest to the median of the final mass distribution. 



fragmented. There are at most one or two density max- 
ima in each one, not many density maxima. These struc- 
tures look much hke t he turbulent cores posited in the 
McKee & Tan (2003) theory for massive star forma- 
tion. When smeared on a resolution of 1700 AU, distinct 
centrally-condensed structures remain visible for three 
of the four massive stars, indicating that these objects 
would be detectable as massive cores in an observation. 
It is important to understand that our analysis says 
nothing about the Lagrangian trajectories of the fluid 
elements that eventually coalesce to form the massive 
stars in our simulations, a topic that has previously re- 



ceiv ed extensive i nvestigat ion by 



and Smith et al. 



Bonnell et al. 



(2004) 



1^' . 

( 2009a|b ), among others. It may well 
be tnat particular fluid elements that are present in the 
cores at the time shown in Figure [TT] do not accrete onto 
the final star and are instead accreted by other stars or 
torn off by turbulent motions, while fluid elements not 
present in the core at the time sh own are eventually ac- 
creted into the final star. Indeed, McKee & Tan (2003) 
predicted in their analytic model that turbulent cores 
should over the course of their lives interact with a sur- 
rounding gas mass comparable to that which eventually 
ends up in their central stars. However, the fact that 
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Fig. 13. — Accretion rate versus stellar mass for the four most 
massive stars present at the end of run TuW. Colored unbroken 
lines indicate the measured simulation accretion rates, while the 
dashed black line is the prediction of the McKee & Tan (2003) 
model (Equation |15| . The simulation accretion rates have been 
smoothed over 50U yr timescales to reduce scatter. 

the Lagrangian elements making up a core change with 
time is irrelevant to the question of whether, as a mas- 
sive star forms, it sits at the center of a gravitationally 
bound Eulerian structure. Figure pT] shows that it does. 
We can make the link between the massive cores and 
the stars they form more quantita tive by comparing to 
the massiv e core evolution niodel of [McKee fc Tan (2002, 
|2003) and |Tan & McKee| (|2004|). I'his model preiJicts 



that the accretion rate onto a star as a function of its 
mass should be 



1.2x10" 



30 Mo 



3/4 



.3/4 / ^* 



■^cl 



1*/ 



1/2 



Moyr- 



(15) 

where m* is the star's instantaneous mass, m*/ is its 
final mass, Sd is the surface density of the molecular 



clump from which it forms, and we have used Mc Kee 
Tanj s fiducial parameter choices, with the exception that 
we have increased the accretion rate by a factor of 2.6 
to in clude subsonic contraction, following [Tan fc McKee| 
([2004). To evaluate this equation and compare it to our 
simulations, we take m^f ~ 10 Mq, since this is roughly 
the mass of our four most massive stars at the end of 
the simulation. For Sd, we note that, in the simulation, 
the core is better-defined than the clump, so we adopt 
McKee fc Ta n's result with Sci replaced by Score- In 
their fiducial model these agree to within a factor :^ 1.2, 
so this does not significantly affect the accretion rate. 
As shown in Figure [TT] our cores have masses of order 
10 Mq in radii of order 0.01 pc, which corresponds to 
Sci = 6.6 g cm~^. With these parameter choices, in 
Figure [13] we plot the accretion rate as a function of 
stellar mass for the four most massive stars at the end 
of the simulation, w hose cores are s hown in Figure [TT] 
and compare to the [McKee fc Tan| prediction. As the 
plot shows, the simulation accretion rates agree quite 
well with the analytic predictions. 

In contrast, the cores that give rise to low mass cores 
(Figure 12) are quite noticeably different from the high 
mass ones. In three of the four cases (the first, third, and 
fourth columns in the Figure) they are also centrally- 
condensed lumps of gas. However, unlike the massive 
cores they are highly sub-fragmented and show many 
density maxima. Clearly these objects are not single 



cores, but instead tightly-packed agglomerations of many 
smaller cores. For the final low mass core (shown in the 
second column of Figure 12) the point at which the star 
forms is a slight overdensity in the middle of a filament, 
and there is no centrally-concentrated object at all. Thus 
massive cores and low mass cores have clearly distinct 
properties. However, we also find that these differences 
are completely indistinguishable in the smeared images, 
indicating that it is not possible to distinguish true high 
mass cores from agglomerations of low mass ones with 
the resolution available in pre- ALMA telescopes, at least 
for objects at the ~kpc distances typical of massive star- 
fo rming refflons . This conclusion is consistent with that 
of lOffner et al.[ ( [2QT21 ). 

It is important to note that the differences between 
high and low mass stars is not simply a function of for- 
mation time. It is certainly true that the most massive 
stars at the end of the simulation preferentially began 
forming early. However, their greater masses are far less 
a refiection of this than it is of their different forma- 
tion environments. The four massive stars grow at time- 
averaged rates of 3.6 — 4.6 x 10~^ Mq yr~^, compared 
to 1.7 - 8.8 X 10"^ Mq yr"^ for the low mass stars. At 
the accretion rates typical of the low mass stars, it would 
require ~ lOtff for one of them to grow to the ~ 10 Mq 
typical of the massive stars. The massive stars are not 
simply those that form first; they are those that form sur- 
rounded by coherent, bound, non-subfragmented struc- 
tures that provide high accretion rates. This is somewhat 
similar to the competitive accretion model in that mas- 
sive stars' preferred locations at the centers of collapsing 
regions that provides their high accretion rates. How- 
ever, it differs from competitive accretion in that these 
cores are non-subgfragmented and have masses at the 
same order of magnitude as the final stars, ^ 10 M©, 
and therefore intermediate between that of the entire 
star cluster, ~ 10^ Mq and the thermal Jeans mass, ~ 1 
Mq. In the competitive accretion model such structures 
should be absent, because everything fragments dow n to 
the thermal Jeans mass (Bonnell et al. 2004; Bate fc Bon-[ 
[iieil||2005| and some objects subsequently grow to larger 
masses by Bondi-Hoyle accretion. There are no ^ 10 
Mq objects that do not subfragment in competitive ac- 
cretion. 

Finally, we note that both the high mass and the low 
mass cores are above the column density threshold S > 1 
g cm~^ for massive star formation posited analytically by 
Krumholz fc McKee ( 2008 ) and confirmed numerically by 
iKrumholz et al., ( ,2010, ) . This means that both the high 
mass and low mass cores have the potential to form mas- 
sive stars; indeed, in one of the four cases shown in Figure 
12 , the low mass star is in fact forming in a core that puts 
most of its mass into a single high mass star. That does 
not appear to be the case for the other three low mass 
stars shown in the Figure, however. Thus a high column 
density is clearly a necessary but not a sufficient condi- 
tion for massive star formation. A high column density 
allows radiative heating to suppress the growth of grav- 
itational instabilities that would lead to fragmentation 
and prevent a massive star from forming. However, if 
the turbulent density field present before a star begins 
radiating is already highly non-linearly fragmented, as 
is the case for several of the low mass cores shown in 
Figure [T2J radiative heating will not undo this fragmen- 
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Fig. 14. — Multiplicity fraction / as a function of system primary 
mass ?77,prim in run TuW. The thick red line shows / as a running 
average. The light red boxes show / computed over discrete bins in 
^^prim- In each case, the width of the box shows the primary mass 
range for that bin, the asterisk shows the mean multiplicity fraction 
for stars in that bin, and the vertical extent of the box shows the 
statistical uncertainty on that value, computed as described in the 
text. Finally, black crosses indicate observational results, with the 
horizontal width indicating the mass range for the observations 
and the vertical range showing the stated uncertainty. The two 
highest mass observational data points are lower limits, indicated 
by the upward arrows. The data shown are taken from, from left 
to right, |Sasri & Reiners (2006) and Allen (2007) (shown as a 
single c omB med pomt), l:^'ischer &; Marc y (1992), Kagh avan et al.| 
([2010), Pr eibisch et al.| ! |iyyyt , and |iVlaso n et al. (2UU9)jthe data 
compilation shown here is the same as that in Bate ( 2012| . 
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Fig. 15. — Semi-major axis versus primary star mass for all the 
binaries in our simulations. For triple and quadruple systems, we 
plot them only once, showing the properties of the most bound pair 
of stars. Points are coded by the mass ratio of the system: purple 
stars for q < 0.1, red squares for g = 0.1 — 0.25, green triangles for 
q = 0.25 — 0.5, and blue circles for q > 0.5. 

tat ion and prevent the core from forming a small cluster 
of low mass stars rather than a few massive ones. 

3.6. Stellar Multiples 

It is also illuminating to consider the properties of the 
stellar multiples that form in run TuW, since producing 
the correct multiplicity fraction has been proposed as 
test for star formatio n models in additio n to producing 
the correct IMF (e.g. [Bonneh et al.|[2QQ7| ). We therefore 
examine the final time slice for this runjj Extracting 

^ In this section of the paper alone we do not exclude stars 
smaller than 0.05 Mq from consideration, but we consider them 
only as companions to larger stars. We allow them to count in this 



the fraction of stars in multiple systems fr om the simula - 
tion requires some care, as pointed out by Bate (2009a). 
Many of the stars in our simulation form a bound clus- 
ter, and thus many stars are bound to many other stars, 
often in hierarchical structures consisting of dozens of in- 
dividual stars; for example a binary and a triple system 
may orbit one another, and these in turn may have ad- 
ditional stars orbiting them. Such agglomerations would 
be extremely unlikely to survive dynamically even for the 
lifetime of a massive star, and would break up if we could 
continue the simulation further. 

Thus we follow Bate in defining stellar multiplicity 
via the following algorithm. We first compute the total 
energy (gravitational plus kinetic in the center of mass 
frame) pairwise for each pair of stars in the simulation. 
We find the most bound system and replace it with a 
single point mass, with a mass equal to the sum of the 
two components, a position located at their center of 
mass, and a momentum equal to the sum of their two 
momenta. We then continually repeat this process, with 
the exception that we do not create aggregates consist- 
ing of more than four individual stars; should the most 
bound system contain five our more stars, we proceed to 
the next most bound pair with fewer than five members 
instead}^ We terminate the process once there are no 
more bound pairs consisting of fewer than five individual 
stars. At the end we are left with a list of star systems, 
some single and some containing up to four individual 
stars. 

Given this list, we can compute the fraction of multiple 
systems as a function of primary star mass. For a set of 
star systems, we define the multiplicity fraction 



f 



B + T + Q 

S + B + T + i 



(16) 



where S^ B, T, and Q and the numbers of single, binary, 
triple, and quadruple systems, respectively^] We choose 
our sets of star systems in two ways. The first is as a 
running average; for a primary mass mprim, we compute 
/ considering all systems for which the primary mass is 
within half a dex of mprim- The second is in discrete 
bins, chosen to roughly match the mass ranges selected 
in observational surveys. We consider primary mass bins 
in the range 0.05 - 0.1 M©, 0.1 - 0.2 M©, 0.2 - 0.5 M©, 
0.5-0.8 Mo, 0.8-1.2 Mo, 1.2-3 Mo, and > 3 Mo- In 
addition to the mean value, we compute the statistical 
uncertainty in this value for each bin 



.put( 





capacity because to omit them would artificially make it impossible 
for stars near our 0.05 Mq cut to have companions. 

^^ The results are not particularly sensitive to the choice of four 
as the maximum size of a system, as long as we stop at some point 
well short of allowing the entire cluster to be considered a single 
large star system. 

11 Following Bat? ('2009a') and 'Rubber fc Whitworth] ( |2005t , 
we measure this quantity rather than either the companion star 
fraction (B + 2T + SQ)/(S + B + T + Q) or the fraction of stars in 
multiple systems {2B + 3T + 4:Q)/(S + 2B + ST + 4Q) because it 
is more robustly determined observationally. If a new member of 
a multiple system is found, for example leading to a binary being 
reclassified as a triple, / does not change, while the companion star 
fraction and the fraction of stars in multiple systems does. 

1^ We determine the statistical uncertainty by assuming that 
there is a true multiplicity fraction /true for stars in that bin, and 
that our sample of systems in that bin represent a series of random 
drawings that follow a binomial distribution. From these assump- 
tions, we can compute the probability distribution for /true given 
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We plot the results in Figure [T4J We see that the 
simulations generally agree well with the observational 
constraints, with the multiplicity fraction reaching near 
unity for stars larger than a few M©, and declining to 
below 0.5 for stars smaller than ~ 0.5 Mq. Our sim- 
ulations somewhat underproduce binaries at the lowest 
masses, which is likely a resolution effect, arising because 
low mass binaries must be very close in order to remain 
bound, and our simulation resolution makes this diffi- 
cult. In our simulation we soften particle-particle grav- 
ity forces on a scale of 0.25 cells, or 5.8 AU, and gas- 
particle gravity forces are necessarily smoothed on the 
grid scale of 23 AU. Thus we cannot easily make binaries 
tighter than ~ 10 AU. At this distance the Keplerian 
speed around a 0.05 Mq object is only 2 km s~^, com- 
parable to the velocity dispersion in the cluster. Thus 
low mass binaries formed in our simulation will tend to 
be broadened and disrupted, and we cannot resolve the 
tighter binaries that will tend to be hardened. This leads 
to an artificial reduction in the bin ary fraction at low 



masses, a phenomenon also noted by [Bate (2012). 

In Figure [15] we illustrate some ot the properties of 
our multiple systems. Systems with more massive pri- 
maries tend to have the smallest separations, as a result 
of dynamical hardening and, in some cases, of a compan- 
ion having been born in the disk of the primary. The 
companions to the most massive stars also tend to be 
fairly massive, with mass ratios of 0.25 — 0.5; this is in- 
consistent with their having been drawn randomly from 
the IMF. These are often triple or quadruple systems. 
Thus we see that the massive stars in our simulation tend 
to form Trapezium- like structures. In contrast, at near- 
Solar masses, the range of semi-major axes and mass 
ratios is extremely broad. 

3.7. Accretion Variability and Outbursts 

A final useful datum to be extract from run TuW is 
the amount of accretion variability for low mass stars, 
which is of interest for its relevance to the protostellar 
luminosity problem and the origin of FU Ori outbursts, 
although due to the duration of our simulation we can 
only address these issues as they apply to class and I 
objects, not class II sources. To characterize the degree 
of luminosity variability we find, we select 12 stars at 
random at the end of our simulation. The final masses 
of these stars range from 0.055 to 1.9 Mq. For each star 
we measure its accretion luminosity, which in our code is 
taken to be 

Lace = 0.75 , (17) 

where m^^ m*, and r* are the star's instantaneous mass, 
accretion rate, and radius, and the factor of 0.75 accounts 
for the energy used to drive protostellar outfio ws (for 
details see the Appendix of Offner et al. 2009). Our 
simulation outputs are spaced roughly 10 — 15 yr apart, 
so this represents the minimum timescale on which we 
can study variability. Outputs are 80 fine grid time steps 
apart, so this timescale is numerically well resolved. 

the measured multiplicity fraction in the simulations /sim and the 
number A^sys of systems in that mass bin. We compute the un- 
certainty by finding the range of values for /true that enclose the 
central 68% of the probability distribution. Note that this range is 
not in general symmetric about /sim- 



In Figure 16 we show the accretion history of each of 
our 12 starsTboth at the maximum temporal resolution 
of the simulations, and smoothed over 200 yr timescales. 
In the figure, the zero of time is the point at which a 
given star forms, and we plot the accretion history for the 
remainder of the simulation. The figure demonstrates 
several interesting results. First, the majority of stars 
have relatively smooth luminosity histories when aver- 
aged over 200 yr timescales. For only a few examples 
do order of magnitude variations in the luminosity occur 
on less than timescales of several kyr. The variability is 
somewhat larger when measured at the maximum tem- 
poral resolution of the simulation, but for most stars this 
is not a large effect. However, there are three exceptions: 
the star shown in blue in the upper right panel, the star 
shown in blue in the lower left panel, and the star shown 
in red in the lower right panel. All of these stars experi- 
ence sudden increases in luminosity on timescales below 
our ability to resolve given the frequency of our output. 
During these spikes, the luminosity rises by 1 — 2 dex 
compared to the long-term average. These are plausibly 
FU Ori- type outbursts, although we caution again that 
these outbursts are occurring in class I sources, not true 
T Tauri stars. 

4. DISCUSSION AND CONCLUSION 

4.1. The Role of Protostellar Outflow Feedback 

In our simulations, we find that protostellar outfiow 
feedback is not particularly effective. Including outfiows 
reduces the star formation rate by only ~ 20%, compara- 
ble to the mass fraction that is ejected from young stars 
by out subgrid model for protostellar winds. This result 
at first a ppears to contradict th ose of previous studies , 

06) 
Matzner (2007), Wang et alj (|j 



at hrst a ppears to contradict tn ose oi previous studies 
including !: [Li fc Nakamur a ( 2006 ) , j Na kamura fc Li|(|2007[) 

''2010 ), and Cunningham 



et al.| ( |2011| ), all of whom find that outfiow lee dback~is 
important. Our results also contradict those of Hansen 
et al.| ( |2012[ ) , who find that outfiow feedback greatly re- 
duces the efficacy of radiative feedback, because it re- 
duces the accretion rate and thus the protostellar lumi- 
nosity. 

Some of our differences from previous results are a 
function of what effects are included in our simulations, 



and in previous work. Wang et al. (2010) note that out- 
flow feedback is only effective as a long-term driver of 
turbulence in the presence of magnetic flelds. Fields fa- 
cilitate transfer of momentum between gas parcels, while 
in their absence most of the momentum injected into a 
protocluster by outflows is simply lost, as the outflows 
break out of the cloud and deposit their momentum and 
energy outside its boundaries. Since our simulations lack 
magnetic flelds, we likely suffer from a similar underesti- 
mate of outffow efficacy. 

A second source of difference is likely to be our choice 
of parameters. We have simulated a fairly massive, high 
surface density cloud representative of the typical Galac- 
tic star-forming region, but with properties quite distinct 
from those of the star- forming regions closest to the Sun. 
All of the previous simulations mentioned above have 
chosen properties typical of these lower density regions. 
Analytic models suggest that protostellar winds are only 
able to eject signiflcant mass from clusters wit h escape 
velocities below i;esc ^ 7 km s~^ (Matzner fc McKee 
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Fig. 16. — Accretion rate versus time for a sample of 12 randomly-selected stars in run TuW. Thick pale lines show accretion luminosities 
averaged over 200 yr timescales, while thin darker lines show the accretion luminosity computed over 10 — 20 yr timescsles, the finest 
available given the frequency with which we output simulation data. There is no distinction between red and blue curves; we simply use 
two different colors to make the two stellar accretion histories shown in each panel more easily distinguishable. Note that the time axes are 
different for the left and right sides. All times are relative to the instant when a star first appears in the simulations, and plots continue 
to the end of the simulation. 



2QQQ[ |Matzner|2QQ7|). In run TuW, the escape velocity is 

i;esc ^ ^GMcjitcl'^) = 4.3 km s"\ within a factor of 2 
of the analytic esti mate. The compar able figure for the 
cluster simulated by [Wang et al.| ( [2QTQ ), which is modeled 
after NG C 1333, i s a tact or ot two lower: ^esc ^2.6 km 



s . For Hansen et al._ (2012), who adopt initial condi- 



tions modeled after p Ophiuchus, it is ^esc ^ 1-6 km s~^. 
Thus it is not surprising that outflows should be much 
more effective in those simulations than in run TuW. In- 
deed, placing the clu ster masses and su rface densities of 
these simulations on Fall et al.| (|2010|)'s diagnostic dia- 
gram for where different sorts of feedback are effective 
(their Figure 2) immediately predicts this dichotomy. 

Given that our simulations likely differ from previous 
work on the importance of outflow feedback due to both a 
physical deficiency (lack of magnetic fields) and a choice 
of parameters that is closer to the typical region than 
most previous work, it is hard to draw general conclu- 
sions about the importance of outfiow feedback in reg- 
ulating star formation. Resolution of this question will 
have to await future magnetohydrodynamic simulations 
that probe the higher density regime we have explored 



in this work. 

4.2. Implications for Massive Star Formation 

The picture of massive star formation that emerges 
from our simulations is generally consist ent with the tur 



bulent core model proposed by McKee & Tan ( 2002 
2003). The massive stars form at the centers of well- 
detined, turbulent, centrally-concentrated structures, 
and these structures feed mass onto the m at a rate that 
is consistent with the predictions of the |McKee fc Tan 
model. In contrast, the regions from which low mass 
stars form are quite noticeably different. They are either 
messy regions consisting of many small density peaks and 
no clear central concentration, or they are small regions 
of filaments. Thus the basic core to star mapping pro- 
posed in the turbulent core model appears to describe 
our simulation fairly well. 

However, we also do see el ements of the altern ative 
competitive accretion model (Bonnell et al.||2007l and 



references therein) operating as well. In particular, our 
massive stars do all form as part of small sub-clusters 
and experience significant dynamical interactions. These 
interactions appear to be important in shaping the multi- 
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plicity properties of the resulting stars, and in producing 
the Trapezium- hke systems in which most of our massive 
stars find themselves at the end of the simulation. 

4.3. Implications for the IMF and the Star Formation 

Rate 

Run TuW represents the first simulation published to 
date that reproduces the observed IMF in a cluster like 
the ONC that is large enough to contain massive stars, 
and where the peak of the mass function is determined 
by a fully self-consistent calculation of gas thermody- 
namics. Previous simulations that have had success in 
reproducing the IMF have either examined small, low- 
density star-forming regio ns that would not be expected 
to produce massive stars (Offner et al. 2009! IBate|2009bl 



2012[ [Urban et al.,2010 . Hansen et al.,2012p , or have re- 
lied on a parameterized, non-selt-consistent equation of 
state to deterniine the location of the IM F peak (e.g. 
Bate fc BonneII||2005[ [Jappsen et al.||2005[ ). 

The success ot run I'uW, in contrast to the failures 
of runs SmNW and TuNW, suggests that obtaining the 
correct IMF from a self-consistent simulation of a typ- 
ical star-forming environ ment is not as simpl e as some 
authors have posited (e.g. Bonnell et al.||2007 ). While a 
lognormal function with a powerlaw tail at high masses 
appears to be a fairly generic result regardless of what 
physics is included in the simulations, getting the peak of 
the lognormal to lie at the correct position in a calcula- 
tion where it is determined self-consistently rather than 
through a hand-imposed equation of state requires care- 
ful attention to the thermodynamics of the gas, which is 
in turn determined primarily by the accretion luminosity 
of protostars. 

This requires several ingredients to work correctly. As 
conjectured in Paper I, the star formation rate cannot 
be too high, and star formation cannot become too cen- 
trally cocentrated; if it is, the resulting accretion lumi- 
nosity becomes so high that formation of low mass stars is 
suppressed and the IMF peak marches to ever-increasing 
mass with time. In addition, outfiows appear to make a 
small but significant contribution by both reducing the 
masses of individual accreting protostars, reducing the 



accretion luminosity, and ejecting mass from the warm 
regions near accreting stars, increasing fragmentation. 
The combination of these effects shifts the mean mass 
downward by a factor of ~ 2. Our results suggest that 
future simulations of gas collapse and fragmentation, if 
they are to reproduce the observed IMF while treating 
the gas thermodynamics self-consistently, must at a min- 
imum include these ingredients. 

We obtain a reduction in the rate and degree of concen- 
tration of star formation in run TuW mainly because we 
have starting our simulations with a self-consistent den- 
sity and velocity field and by embedding the protocluster 
gas clump in a realistic surrounding turbulent molecular 
cloud that continually feeds energy into it via a turbu- 
lent cascade, rather than simulating a smooth, isolated 
clump as in most previous work. However, the star for- 
mation rate per free-fall time we obtain is still an order of 
magnitude too high compared to observations, likely as 
a result of other mechanisms we have not included. For 



examp le, Price & Bate (2009) and Padoan & Nordlund 
(2011 ) both find that magnetic fields with strengths com- 
parable to those in observed molecular clouds reduce star 
formation rates by a factor of a few to ^ 10. Depend- 
ing on protocluster properties, ionized gas and radiation 
pressure may also contribute to reducing the star forma- 
tion rate (e.g. [Fah et"aL]|2010| [Lopez et aLpOllj ). Ulti- 
mately, since accretion luminosity plays a critical role in 
regulating the IMF, the problems of determining the star 
formation rate and the IMF cannot be fully separated. 
A truly accurate simulation must reproduce both. 
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